module params_mod
  use spark_cubed_sphere, only: proc, r_kind, rad, D2R
  use flogger
  use datetime
  use time_mod, only: start_time, end_time

  implicit none

  integer         :: case_id
  character(99)   :: case_name                  = 'hope_output'
  integer         :: run_days
  integer         :: run_hours
  integer         :: run_minutes
  integer         :: run_seconds
  real(r_kind)    :: dt
  real(r_kind)    :: history_interval           = 86400
  real(r_kind)    :: print_interval             = 3600
  integer         :: hdim                       = 2
  integer         :: vdim                       = 3
  integer         :: nx                         = 0
  integer         :: ny                         = 0
  integer         :: nz                         = 1
  integer         :: nc                         = 0
  real(r_kind)    :: dx, dy, deta
  integer         :: model_type                 = 1         ! 1 for shallow water, 2 for baroclinic core
  character(99)   :: time_scheme                = 'SSP_RK3' ! Choose from 'Euler', 'SSP_RK3', 'RK4'
  character(99)   :: riemann_scheme_hor         = 'AUSM_up' ! Horizontal Riemann Solver, choose from 'AUSM_up', 'LLF'
  character(99)   :: riemann_scheme_ver         = 'AUSM_up' ! Vertical Riemann Solver, choose from 'AUSM_up', 'LLF'
  character(99)   :: recon_h_scheme             = 'WENO2D'  ! Choose from 'WENO1D', 'WENO2D', 'WENO2D_rhombus', 'Poly1D', 'Poly2D', 'Poly2D_rhombus'
  integer         :: recon_h_order              = 5
  character(99)   :: recon_v_scheme             = 'Poly1D'
  integer         :: recon_v_order              = 5
  logical         :: use_tci                    = .true.
  integer         :: IRK_stage                  = 5
  real(r_kind)    :: IRK_residual               = 1.E-5
  integer         :: IRK_max_iter               = 30
  real(r_kind)    :: damp_thickness             = 10000
  real(r_kind)    :: damp_max_coef              = 1
  
  ! Baroclinic only
  character(100) :: planet = 'earth'
  character(100) :: topo_type
  character(999) :: topo_file
  integer        :: v_dist_type = 1
  integer        :: v_coord_type= 2
  real(r_kind)   :: m_coef = 2.5       ! Vertical coordinate streching coef
  real(r_kind)   :: max_height = 40000 ! Maximum vertical height
  logical        :: use_3d_gravity = .false.
  
  ! Derivative parameters
  integer nsubstep
  integer sw
  integer swv
  integer hw
  integer rw
  integer rwv
  integer neq
  integer nfq
  integer nfqv
  integer ncq
  integer nvq
 
  integer :: nVar = 3
  integer :: ndom = 6
  integer :: dvs, dve ! Rayleigh damping variable start/end indices

  namelist /hope_control/    &
    case_id,                 &
    case_name,               &
    start_time,              &
    end_time,                &
    run_days,                &
    run_hours,               &
    run_minutes,             &
    run_seconds,             &
    dt,                      &
    time_scheme,             &
    riemann_scheme_hor,      &
    riemann_scheme_ver,      &
    history_interval,        &
    print_interval,          &
    hdim,                    &
    vdim,                    &
    ndom,                    &
    nc,                      &
    nx,                      &
    ny,                      &
    nz,                      &
    model_type,              &
    recon_h_scheme,          &
    recon_h_order,           &
    recon_v_scheme,          &
    recon_v_order,           &
    use_tci,                 &
    IRK_stage,               &
    IRK_residual,            &
    IRK_max_iter,            &
    damp_thickness,          &
    damp_max_coef,           &
    ! For barolinic initialize
    planet,                  &
    topo_type,               &
    topo_file,               &
    v_dist_type,             &
    v_coord_type,            &
    m_coef,                  &
    max_height,              &
    use_3d_gravity

contains

  subroutine params_init(file_path)

    character(*), intent(in) :: file_path
    
    type(timedelta_type) delta_t
    
    open(10, file=file_path)
    read(10, nml=hope_control)
    close(10)
    
    if (nc /= 0) then
      nx = nc; ny = nc
    else if (nx /= 0 .and. ny /= 0) then
      if (ny /= nx) then
        call log_warning('nx /= ny, dy is reset to nx!', pid=proc%id)
        ny = nx
      end if
      nc = nx
    end if

    dx = 90.d0 / nx
    dy = dx
    
    dx = dx * D2R
    dy = dy * D2R
    deta = max_height / nz
 
    select case (trim(time_scheme))
    case ('Euler')
      nsubstep = 2
    case ('SSP_RK3')
      nsubstep = 3
    case ('WRF_RK3')
      nsubstep = 3
    case ('RK4')
      nsubstep = 4
    case ('IMEX')
      nsubstep = 2
    end select
    
    if(model_type==0)then
      nVar = 1
      nz = 1
      hdim = 2
      vdim = 1
      swv = 0
      rwv = 0
    elseif(model_type==1)then
      nVar = 3
      nz = 1
      hdim = 2
      vdim = 1
      swv = 0
      rwv = 0
    elseif(model_type==2)then
      nVar = 6
      !hdim = 3
      if(hdim==3)then
        call log_notice('hdim==3, reset recon_v_order to 3', pid=proc%id)
        recon_v_order = recon_h_order
      endif
      swv = recon_v_order
      rwv = ( recon_v_order - 1 ) / 2
      if( trim(time_scheme) == 'IMEX' )then
        call log_notice('vdim is reset to 1 for IMEX', pid=proc%id)
        vdim = 1
      endif
    endif
   
    sw = recon_h_order
    rw = (sw - 1) / 2
    hw = rw + 1
    select case (recon_h_scheme)
    case ('WENO1D', 'Poly1D')
      neq = 1
    case ('WENO2D', 'WENO2D_rhombus', 'Poly2D', 'Poly2D_rhombus', 'Poly3D')
      neq = ( recon_h_order + 1 ) / 2
    end select
    
    if(hdim==2)then
      ncq  = neq**2
      nfq  = neq
      nvq  = ncq
    elseif(hdim==3)then
      ncq  = neq**2
      nfq  = neq**2
      nvq  = neq**3
    endif
    
    if(vdim==1)then
      nfqv = 1
    elseif(vdim==3)then
      nfqv = ncq
    else
      print*,'Choose vdim from 1 or 3, vdim=',vdim
      stop
    endif
    
    if( .not. ( trim(recon_h_scheme)=='WENO2D' .or. trim(recon_h_scheme)=='WENO2D_rhombus' ) )use_tci = .false.

    if (end_time < start_time) end_time = start_time
    if (run_days    > 0)then
       call delta_t%init(days   =run_days   )
       end_time = end_time + delta_t
    endif
    if (run_hours   > 0)then
      call delta_t%init(hours  =run_hours  )
      end_time = end_time + delta_t
    endif
    
    if (run_minutes > 0)then
      call delta_t%init(minutes=run_minutes)
      end_time = end_time + delta_t
    endif

    if (run_seconds > 0)then
      call delta_t%init(seconds=run_seconds)
      end_time = end_time + delta_t
    endif
    
  end subroutine params_init

end module params_mod
